# 相关库
from scipy import  stats
from statsmodels.graphics.tsaplots import plot_acf as ACF   # 自相关图
from statsmodels.graphics.tsaplots import plot_pacf as PACF   # 偏自相关图
from statsmodels.tsa.stattools import adfuller as ADF # 平稳性检测 
from statsmodels.tsa.ar_model import AutoReg
from statsmodels.tsa.arima.model import ARIMA
from statsmodels.graphics.api import qqplot
from scipy.stats import shapiro
import statsmodels.api as sm  # 统计相关的库
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import arch  # 条件异方差模型相关的库
import seaborn as sns     # seaborn画图
sns.set(color_codes=True) # seaborn设置背景

#①导入数据
data=pd.read_excel('物联网板块价格截取.xlsx')
data.set_index('date', inplace=True) #设定日期为索引

#r2=np.log(data['close'])-np.log(data['close'].shift(1)) #计算对数收益率
#r2=r2.dropna()  #去除包含NaN的行

d1 = data['close']
r1 = np.array(d1)
t = ADF(r1)
output=pd.DataFrame(index=['Test Statistic Value', "p-value", "Lags Used", "Number of Observations Used","Critical Value(1%)","Critical Value(5%)","Critical Value(10%)"],columns=['value'])
output['value']['Test Statistic Value'] = t[0]
output['value']['p-value'] = t[1]
output['value']['Lags Used'] = t[2]
output['value']['Number of Observations Used'] = t[3]
output['value']['Critical Value(1%)'] = t[4]['1%']
output['value']['Critical Value(5%)'] = t[4]['5%']
output['value']['Critical Value(10%)'] = t[4]['10%']
print(output)

print("收盘价序列ADF检验平稳性结果：")
print("p-value:   ",t[1])
print("\n")

d2=data['close'].diff()
r2=np.array(d2)[1:]
t = ADF(r2)
print("收盘价一阶差分序列ADF检验平稳性结果：")
print("p-value:   ",t[1])
print("\n")

d1.plot(figsize=(15,4)) #绘制收盘价序列折线图
print("绘制收盘价序列折线图")
print("\n")
plt.show()

#②使用ACF和PACF对AR和MA模型定阶
#ACF定阶
fig = ACF(r1, lags = 30)
print("绘制ACF图")
print("\n")
plt.show()
#PACF定阶
fig = PACF(r1, lags = 30)
print("绘制PACF图")
print("\n")
plt.show()

#信息准则定阶（运行速度慢）
#print(sm.tsa.arma_order_select_ic(r1,max_ar=10,max_ma=10,ic='aic')['aic_min_order'])  # AIC

#③ARIMA建模
temp = np.array(r1) # 载入收盘价序列
model =ARIMA(temp,order=(1, 0, 8))  
res = model.fit()  
print(res.summary())

plt.rcParams['font.sans-serif'] = ['simhei'] #字体为黑体
plt.rcParams['axes.unicode_minus'] = False #正常显示负号 
plt.figure(figsize=(10,4))
plt.plot(temp,'b',label='收盘价序列')
plt.plot(res.fittedvalues, 'r',label='ARIMA model')
plt.legend()
print("绘制ARIMA模型")
print("\n")
plt.show()

#④残差白噪声检验
delta = res.fittedvalues  - temp  # 残差
plt.figure(figsize=(10,6))
plt.plot(delta,'r',label=' residual error')
plt.legend(loc=0)
acf,q,p = sm.tsa.acf(delta,nlags=12,qstat=True)  ## 计算自相关系数 及p-value
out = np.c_[range(1,13), acf[1:], q, p]
output=pd.DataFrame(out, columns=['lag', "AC", "Q", "P-value"])
output = output.set_index('lag')
print("残差白噪声检验结果：")
print("\n")
print(output)

#⑤拟合优度
score = 1 - delta.var()/temp.var() #计算拟合优度
print("拟合优度：")
print(score)

# 指定预测步数
steps = 5  # 假设预测未来5个时间点

# 使用已拟合的ARIMA模型进行预测
forecast = res.forecast(steps=steps)

# 获取预测结果
forecast_series = pd.Series(forecast, index=range(len(temp), len(temp) + steps))
print("未来{}个时间点的收盘价预测：".format(steps))
print(forecast)

# 可视化预测结果
plt.figure(figsize=(10, 6))
plt.plot(temp, label='历史收盘价')
plt.plot(forecast_series, color='red', linestyle='--', label='预测收盘价')
plt.xlabel('时间')
plt.ylabel('收盘价')
plt.title('收盘价预测')
plt.legend()
plt.show()


